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Abstract 

This  paper  begins  with  an  overview  of  the  boundary  condition  cap¬ 
turing  approach  to  solving  problems  with  interfaces.  Although,  the 
authors’  original  motivation  was  to  extend  the  ghost  fluid  method 
from  compressible  to  incompressible  flow,  the  elliptic  nature  of  in¬ 
compressible  flow  quickly  quenched  the  idea  that  ghost  cells  could 
be  defined  and  used  in  the  usual  manner.  Instead  the  boundary 
conditions  had  to  be  implicitly  captured  by  the  matrix  formulation 
itself,  leading  to  the  novel  approach.  We  first  review  the  work  on 
the  variable  coefficient  Poisson  equation,  noting  that  the  simplic¬ 
ity  of  the  method  allowed  for  an  elegant  convergence  proof.  Sim¬ 
plicity  and  robustness  also  allowed  for  a  quick  extension  to  three- 
dimensional  two-phase  incompressible  flows  including  the  effects  of 
viscosity  and  surface  tension,  which  is  discussed  subsequently.  The 
method  has  enjoyed  popularity  in  both  computational  physics  and 
computer  graphics,  and  we  show  some  comparisons  with  the  tradi¬ 
tional  delta  function  approach  for  the  visual  simulation  of  bubbles. 
Finally,  we  discuss  extensions  to  problems  where  the  velocity  is  dis¬ 
continuous  as  well,  as  is  the  case  for  premixed  flames,  and  show 
an  example  of  multiple  interacting  liquids  that  includes  all  of  the 
aforementioned  phenomena. 


1  Introduction 

Enforcing  a  variety  of  boundary  conditions  or  jump  conditions  at  interfaces 
is  important  for  developing  accurate  numerical  methods  in  the  held  of  com¬ 
putational  physics.  For  example,  Fedkiw  et  al.  [1]  developed  the  ghost  fluid 
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method  to  capture  the  boundary  conditions  at  a  contact  discontinuity  in 
the  inviscid  Euler  equations.  This  method  was  later  extended  to  treat  more 
general  compressible  flow  discontinuities  such  as  shocks,  detonations,  and 
deflagrations  in  [2].  For  a  fully  conservative  implementation,  see  [3]. 

Motivated  by  the  ghost  fluid  method,  [4]  developed  a  novel  boundary 
condition  capturing  approach  for  solving  the  variable  coefficient  Poisson 
equation  in  the  presence  of  interfaces  where  both  the  variable  coefficients 
and  the  solution  itself  may  be  discontinuous.  The  method  is  robust  and 
easy  to  implement  even  in  three  spatial  dimensions.  Furthermore,  the  co¬ 
efficient  matrix  of  the  associated  linear  system  is  the  standard  symmetric 
matrix  for  the  variable  coefficient  Poisson  equation  in  the  absence  of  in¬ 
terfaces  allowing  for  straightforward  application  of  standard  “black  box” 
solvers.  A  formal  convergence  proof  was  given  in  [5]  based  on  the  weak 
formulation.  Both  the  method’s  simplicity  and  its  three-dimensional  for¬ 
mulation  made  it  an  ideal  candidate  for  three-dimensional  flows,  and  [6] 
extended  the  method  to  treat  two-phase  incompressible  flow  including  the 
effects  of  viscosity,  surface  tension  and  gravity.  While  the  traditional  finite 
difference  techniques  for  two-phase  incompressible  flow  involved  numeri¬ 
cal  smearing  of  the  discontinuous  quantities  near  the  interface,  see  e.g. 
[7,  8,  9],  [6]  treats  the  interface  in  a  sharp  fashion  dramatically  reducing 
parasitic  currents.  Interestingly,  [10]  has  shown  that  these  delta  function 
smearing  techniques  may  lead  to  0(1)  errors  in  even  basic  computations 
if  used  in  the  standard  fashion  (see  also  [11,  12]).  While  there  have  been 
other  sharp  interface  methods  such  as  the  immersed  interface  method  [13], 
complicated  implementations  and  issues  with  stability  and  accuracy  have 
precluded  their  use  in  some  of  the  more  complex  problems.  Notably,  [14] 
proposed  a  new  simpler  form  of  the  immersed  interface  method  that  is  quite 
similar  to  the  boundary  condition  capturing  approach  of  [4]. 

The  methods  proposed  in  [4]  and  [6]  have  also  enjoyed  popularity  out¬ 
side  the  computational  physics  community.  For  example,  [15]  indepen¬ 
dently  used  these  techniques  for  the  visual  simulation  of  incompressible 
viscous  two-phase  fluids  with  realistic  small-scale  details.  Although  they 
used  a  slightly  simplified  approach,  decoupling  the  pressure  discontinuity 
from  the  viscosity  discontinuity,  they  still  obtained  sharp  interface  profiles 
of  solution  variables  and  showed  that  the  method  is  highly  preferable  to  a 
smeared  out  delta  function  approach.  Similar  techniques  were  also  used  to 
produce  realistic  three-dimensional  simulations  of  fire  in  [16].  This  work 
was  based  on  the  computational  physics  research  of  [17]  that  extended  the 
methods  proposed  in  [4]  and  [6]  to  treat  two-phase  incompressible  flow 
problems  where  one  phase  is  being  converted  into  another,  e.g.  the  burn¬ 
ing  of  a  premixed  flame  or  the  vaporization  of  liquid  water.  Most  recently, 
[18]  combined  all  of  these  techniques  into  a  single  simulation  framework  for 
multiple  interacting  liquids. 
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2  Variable  Coefficient  Poisson  Equation 

Based  on  the  ghost  fluid  method,  [4]  used  a  boundary  condition  capturing 
approach  to  develop  a  new  numerical  method  for  the  variable  coefficient 
Poisson  equation  in  the  presence  of  interfaces  where  both  the  variable  co¬ 
efficients  and  the  solution  itself  may  be  discontinuous.  The  method  gives  a 
simple  dimension  by  dimension  discretization  that  can  be  readily  applied 
in  three  spatial  dimensions.  The  resulting  linear  system  is  the  standard 
symmetric  discretization  for  the  Poisson  equation,  allowing  for  fast  solvers 
to  be  used.  Notably,  the  method  maintains  a  sharp  profile  and  does  not 
suffer  from  numerical  smearing  at  the  interface  unlike  approaches  based  on 
a  delta  function  formulation. 

Consider  a  Cartesian  computational  domain,  Q,  with  exterior  boundary, 
dCl,  and  a  lower  dimensional  interface,  T,  that  is  defined  by  a  level  set 
function,  <j>,  and  divides  the  computational  domain  into  disjoint  pieces,  f l~ 
and  fi+.  The  variable  coefficient  Poisson  equation  is  given  by 

V  •  (0(x)Vu(x))  =  /(£),  f  eO  m 

u(x)  =  g(x),  x  €  dfl 


where  x  =  ( x,y,z )  are  the  spatial  dimensions,  V  =  (gy>  gy)  is  the 
gradient  operator,  and  j3{x)  is  presumed  to  be  continuous  on  each  disjoint 
subdomain,  and  f 1+ .  The  jump  conditions  or  internal  boundary  condi¬ 
tions  are  specified  along  the  interface  T  as 


where 


Mr  =  a(x), 

x  e  r 

(2) 

[ Kir  =  b(%)> 

x  g  r 

Mr  =  u+(x)  — 

u~(x) 

(3) 

[Pun\ r  =  (3+(x)u+(x) 

-  /3-(x)u~{x) 

specifies  the  direction  of  the  jump  with  the  “±”  subscripts  referring  to  fi  , 
Here  un  =  X7u-N  is  the  normal  derivative  of  u  with  N  the  local  unit  normal 
to  the  interface.  The  use  of  Dirichlet  boundary  conditions  in  equation  1  is 
for  exposition  only  and  Neumann  boundary  conditions  of  un(x)  =  g{x)  for 
x  €  dfi  could  be  used  on  the  outer  boundary  instead. 


2.1  Incorporating  Jump  Conditions 

For  simplicity,  consider  the  one-dimensional  variable  coefficient  Poisson 
equation 


((3ux)x  =  f(x) 


(4) 
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with  fixed  Dirichlet  boundary  conditions  on  dfl  and  a  standard  second 
order  discretization  of 


Uj  +  1—Uj  \  _  a  (  Uj-Uj-t  \ 

Ai  J  Pi- 1  ^  Ax  j 


Ax 


=  fi 


(5) 


for  each  unknown  ut .  At  the  fluxes,  /3i±i  =  /3(xi±i)  are  defined  in  accor¬ 
dance  with  the  side  of  the  interface  the  flux  is  located  on  as  determined  by 


If  Xk  €  and  Xk+i  €  fi+,  one  can  see  that  Uk+^xUk  is  O(^),  while 
all  the  other  terms  of  the  form  u%+^  u%  are  0(1)  and  approximate  the 


local  derivative.  The  term  Ufc+.!  Uk  must  be  modified  in  accordance  with 

Ate 

the  jump  condition  in  order  to  obtain  a  reasonable  approximation  of  the 
derivative  near  the  interface.  In  particular,  one  should  avoid  naively  mixing 
terms  from  different  domains.  Drawing  inspiration  from  the  ghost  fluid 
method,  one  can  define  wjT+1  =  ttfc+i  —  ar  and  u ^  =  uk  +  ar,  where 


_  dfcl^fc+il  +  afc+i|^fci 

or  “  Wt\  +  ¥kA\  (6) 

gives  the  jump  at  the  interface.  Then  in  the  equation  for  grid  node  k,  we 
replace  uk+i  with  uAv  to  obtain 


(uk+i  —  ar)  —  uk  \  o  ( Uk—Uk- 1  ' 

Ax  J  Pk- 1  ^  Ax  J 


Ax 


fk- 


(7) 


Similarly,  the  equation  for  grid  node  k  +  1  is  modified  to  be 

pk+i  (Mfc+2AT+1)  -4+j 


Ax 


Ax 


f k-\- 1  • 


(8) 


Next,  we  turn  our  attention  to  the  derivative  jump  condition  where  the 
precise  location  of  the  interface  plays  a  larger  role,  and  thus  we  use 


hftfel 

I'/’fcl  +  l^fc+il 


(9) 


to  estimate  the  subcell  interface  location.  That  is,  the  interface  splits 
the  cell  into  two  pieces  of  size  9 Ax  on  the  left  and  (1  —  6)  Ax  on  the 
right.  Denoting  the  value  of  u  at  this  subcell  interface  location  by  ui,  and 
interpolating  the  derivative  jump  at  the  interface  as 

,  bk\(j)k+1\  +  bk+i\<t>k\  /1flx 

br  -  M  +  fc+.l  '  (10) 
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one  can  discretize  the  derivative  jump  condition  as 


q+  /  uk+ 1  -  ui  \ 
and  solve  for  uj  as 


/ 3+uk+\9  +  0  uk(l  —  9)  —  brO(l  —  9) Ax 
0+0  +  0- (1-0) 


(11) 


(12) 


so  that  approximations  to  the  derivatives  on  the  left  and  right  sides  of  the 
interface  can  be  written  as 


/ -  uk\  _  ~  / uk+ 1  -uk\  0br(l  -  9) 

^  oax  y  ;  V  A®  )  P+ 


+  (  Uk+ 1  -Ul  \  _  a  f  Uk+l  -Uk\  0br9 
/  \  (1  —  9) Ax )  V  A^  )  + 

where 

0  = _ _ 

^  0+6  +  0-(l-9) 

defines  an  effective  0.  That  is,  we  have 


K 


(•^fc+i  —  Qr)~  Uk 

Ax 


br(l-e) 

0+ 


Ax 


=  fk 


and 


Uk  +  2  —  Uk+1 
Ax 


Uk+i  —  juk+ar) 
Ax 


Ax 


fk+1 


(13) 

(14) 

(15) 

(16) 

(17) 


as  the  equations  for  the  unknowns  uk  and  uk+ i  respectively. 
Of  course,  these  cau  be  rewritten  as 


a  (  Wfc+i—  Uk\  a  (  Uk—Uk- i  \ 

^  (v — sit — ;  -  /5fc-i  {  Ax  ) 


Ax 


=  fk 


0ar  0br(l  -  9) 
(Ax)2  +  0+  Ax 


(18) 


and 


O  (  Uk  +  2  —  Uk+l\  O  ( Uk+1—Uk\ 

Pk+  §  ^  &j  )  P  \  Ax  ) 


Ax 


—  fk+1  — 


0ar 

(Ax)2 


f3br9 
0~  Ax 


(19) 
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Figure  1:  {(3ux)x  =  f(x),  [u]  ±  0,  \(3un]  ^  0 


to  emphasize  that  this  numerical  method  yields  a  symmetric  linear  system 
with  (3k+i  =  (3. 

In  figure  1,  (3  =  2  and  f(x)  =  (8x2  —  4)e~x  on  the  interior  region, 
and  (3=1  and  f(x)  =  0  on  the  exterior.  At  x  =  .3,  [u]  =  — e-09  and 
[flun]  =  —  1.2e-09,  while  at  x  =  .6,  [u]  =  — e-  36  and  [(3un\  =  2.4e-  36. 
The  solution  is  plotted  on  top  of  the  exact  solution  of  u(x)  =  e~x  on 
the  interior  and  u(x)  =  0  on  the  exterior.  Note  that  the  sharp  jumps  are 
preserved  in  both  the  function  and  its  derivatives. 

2.2  Multiple  Spatial  Dimensions 

Consider  the  two-dimensional  Poisson  equation 

(f3ux)x  +  ( (3uy)y  =  f{x)  (20) 

with  interface  jump  conditions,  [u]r  =  a(xp)  and  [f3un] y  =  6(xr)-  The  unit 
normal  is  N  =  ( n1,n 2)  with  (f>  <  0  in  and  <f>  >  0  in  Q+  implying  that 
the  unit  normal  points  from  f into  f 1+. 

The  normal  and  tangential  derivatives  can  be  defined  in  terms  of  ux, 
uy  and  N  as 

un  =  uxrr  +  uyn 2  (21) 
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and 


Ut  =  uxn2  —  uyii 1  (22) 

respectively.  Then 

ux  —  unn 1  +  Utn2  (23) 

and 

=  unn 2  —  Utn1  (24) 

follow  directly  from  equations  21  and  22.  Multiplying  equations  21  and  22 
by  (3  and  taking  the  jump  across  the  interface  leads  to 

\(3un\  r  =  {Pu^rn1  +  [ Puv\Tn 2  (25) 

and 

[/3'«t]r  =  [/3'Ua-lrn2  -  [pUy^n1  (26) 

noting  that  N  is  continuous  across  the  interface.  In  the  same  fashion, 

\fiux\r  =  [ dun\Tnl  +  [P'ut}rn2  (27) 

and 

[Puv] r  =  [f3un]rn2  -  [pu^rn1  (28) 

can  be  obtained  from  equations  23  and  24. 

Suppose  that 

\Pux\t  =  [P'Unjrn1  (29) 

and 

[Puy\ r  =  [ Pun\rn 2  (30) 

are  used  in  place  of  equations  27  and  28.  While  equations  29  and  30  are 
false  in  general,  they  still  lead  to  an  identity  when  plugged  into  equation 
25.  However,  they  lead  to  [Pup r  =  0  when  plugged  into  equation  26. 
That  is,  equations  29  and  30  allow  one  to  correctly  capture  the  jump  in 
the  normal  derivative  while  smearing  out  the  jump  in  the  tangential  deriv¬ 
ative.  More  importantly,  equations  29  and  30  allow  the  derivative  jump 
condition,  [Pun] r  =  fr(xr),  to  be  rewritten  as  two  separate  jump  condi¬ 
tions,  [ Pux\t  =  b(x rjn1  and  [Puy\ r  =  b(xr)n2,  allowing  a  dimension  by 
dimension  application  of  the  numerical  method. 
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Figure  2:  V  •  (f3Vu)  =  f(x,y ),  [u]  fy  0,  [/3un]  fy  0 


In  two  dimensions,  the  equation  is  discretized  at  each  grid  point  (*,  j) 
as 


'"'a/"')  fy,  ■(""  V  ) 


A  y 


fi,j  +  Fx  +  Fy 


(31) 


and  included  in  the  linear  system  of  equations.  Here  Fx  =  FL  +  FR  and 
Fy  =  Fb+Ft  ,  where  FL ,  FR,  FB  and  FT  incorporate  the  jump  conditions 
across  the  left,  right,  bottom  and  top  arms  of  the  stencil  respectively. 

Figure  2  shows  an  interface  is  defined  by  ( x(9),y(9 ))  where  x(6)  = 
.6cos(0)  —  .3cos(30),  y(9)  =  1.5  +  .7sin(0)  —  .O7sin(30)  +  .2sin(70),  6  S 
[0,  27t).  The  exact  solution  is  u(x,y)  =  ex(x2sm(y)  +  y2)  on  the  interior 
and  u(x,  y)  =  — ( x 2  +  y2)  on  the  exterior.  (3  =  1  with  f(x,  y)  =  ex(2  +  y2  + 
2sin(y)  +  4xsin(y))  on  the  interior  and  (3  =  10  with  f(x,y)  =  —40  on  the 
exterior.  The  jump  conditions  are  [u]  =  — ( x 2  +  y2)  —  ex(x2  sin(y)  +  y2) 
and  [ (3un ]  =  (—20a’  —  ex((x2  +  2x)  sin (y)  +  y2))n i  +  (— 20y  —  ex(x2  cos (y)  + 
2y))ri2-  Note  that  the  sharp  jumps  are  preserved  in  both  the  function  and 
its  derivatives  even  in  multiple  spatial  dimensions. 

The  three-dimensional  scheme  follows  a  similar  dimension  by  dimension 
framework  with  \(3ux\t  =  [/3'Un]r?i1,  [/ 3uv]r  =  [/3un\rn2,  and  \(3uz]r  = 
[/3zi„]rn3  false,  but  still  correctly  specifying  the  normal  derivative  jump 
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condition  while  smearing  out  the  tangential  derivatives. 


3  Multiphase  Incompressible  Flow 

[6]  adapted  this  boundary  condition  capturing  technique  to  three-dimensional 
multiphase  incompressible  flow  calculations  including  the  effects  of  viscos¬ 
ity,  surface  tension  and  gravity.  They  used  a  projection  method  to  make  an 
intermediate  velocity  field  divergence  free  at  each  time  step,  thus  requiring 
the  solution  of  a  variable  coefficient  Poisson  equation  with  a  discontinuous 
solution. 

The  basic  equations  for  viscous  incompressible  flow  are 

pt  +  V  •  V/9  =  0  (32) 


(v-r)T 

P 


+  9 


(33) 


where  t  is  the  time,  p  is  the  density,  V  —<  u,v,w  >  is  the  velocity,  p  is  the 
pressure,  p  is  the  viscosity,  g  =<  0,g,  0  >  is  gravity,  and  r  is  the  viscous 
stress  tensor, 


T=  P 


2  ux 

Uy  +  VX 

Uz  +  wx  \ 

\  (  Vu 

\  f 

Uy  +  VX 

2  vy 

VZ  +  Wy 

J  =  p  j  Vw 

+  p\  Vv 

uz  +  wx 

VZ  +  Wy 

2  wz  ) 

1  V  Vw  ) 

1  \  X7w 

.(34) 


These  equations  are  trivially  derived  from  the  viscous  compressible  Navier- 
Stokes  equations  using  the  divergence  free  condition,  V  •  V  =  0. 

Following  the  projection  method,  an  intermediate  velocity,  V* ,  is  cal¬ 
culated  via 


and  then  the  velocity  field  at  the  new  time  step  is  obtained  through 


vn+ 1  -  V *  v» 

- 7 - h  —  =  0. 

At  p 

Taking  the  divergence  of  equation  36  results  in 


V  • 


V  •  V * 
At 


(36) 


(37) 


after  setting  V  •  Vn+1  to  zero. 
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3.1  Jump  Conditions 

The  boundary  condition  capturing  approach,  and  sharp  interface  approaches 
in  general,  require  a  set  of  jump  conditions.  Although  nontrivial  for  the 
three-dimensional  Navier-Stokes  equations  including  the  effects  of  viscos¬ 
ity,  surface  tension,  and  gravity,  they  were  derived  in  [6].  We  repeat  that 
derivation  here. 

Applying  conservation  allows  one  to  write  the  jump  conditions  for  an 
interface  moving  with  the  local  fluid  velocity  in  the  normal  direction  as 

^  fy  j  ( pI-T)NT 


where  Tj  and  T2  are  orthogonal  unit  tangent  vectors,  I  is  the  identity 
matrix,  a  is  the  coefficient  of  surface  tension  (a  constant),  and  k  =  — V  •  N 
is  the  local  curvature  of  the  interface.  Equation  38  states  that  the  net  stress 
on  the  interface  must  be  zero  (since  it  has  no  mass).  For  more  details,  see 
[19,  20], 

Using  the  definition  of  r  in  equation  38  leads  to 


/  p\  (  N  \  /  Vic  N  \ 
0  -  /x  fi  Vi>  •  N 

\  0  /  \  T2  J  \  Vw  ■  N  J 


/  Vu 

■N 

Vu 

•  N 

Vw 

■N  \ 

" 

(  OK  \ 

Vjc 

T\ 

Vv  • 

Ti 

Vru  ' 

■Ti 

•  N 

=  0 

V  V  vi  ■ 

Ti 

Vu  ■ 

T2 

Vw  ' 

Ti  / 

V  0  I 

which  can  be  written  as  three  separate  jump  conditions 


2/x( 


p  —  2/i  Vw  •  A/",  \7v  •  N,  \7w  •  N )  •  N 


=  an 


(40) 


fi  ( Vw  •  N )  Vv  •  N,  \7w  •  N)  •  Ti  + 


fi  (  Vu  •  Ti,  \7v  •  Ti,  \7w  •  T\ )  •  N 


[i  ( Vu  •  N ,  Vv  •  N,  Vw  •  N )  •  T2+ 


fi  [Vu  •  T2,  Vv  •  T2,  \7w  •  r2J  •  N 
Since  the  flow  is  viscous,  the  velocities  are  continuous 
\u]  =  \v]  =  \w]  =  0 


=  0  (41) 

=  0  (42) 

(43) 
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as  well  as  their  tangential  derivatives 


[Vu  •  T{\  =  [Vv  ■  fi]  =  [V«j  •  fi]  =  0 


(44) 


[Vu  •  f2]  =  [Vu  •  f2]  =  [Vw  •  f2]  =  0  (45) 

so  that  the  identity 

(Vu  •  N,  Vv  ■  N,  Vw  •  N^j  ■  N  +  (Vu  ■  Ti,  Vv  ■  fu  Vw  •  7\)  •  Tx+ 

•  f2,  Vu  ■  T2,  Vw  ■  f2^j  ■  T2  =  V  •  V  =  (j)46) 

can  be  used  to  obtain 

(Vu  •  N,  Vv  ■  N,  Vw  •  ■  N  =0  (47) 

emphasizing  that  the  normal  derivative  of  the  normal  component  of  the  ve¬ 
locity  is  continuous  across  the  interface  allowing  equation  40  to  be  rewritten 
as 

[p]  -  2  [p]  (Vu  ■  N,Vv  ■  N,  Vw  ■  Jv'j  ■  N  =  <jk  (48) 

Next,  the  family  of  identities  of  the  form 

[AB\  =  B[A]  +  A[B]  (49) 


A  —  ciArjjgjli  -\-  bAiejti  B  —  bBright  -t-  ciBieft,  a  -\-  b  —  1  (50) 

is  used  along  with  equations  44  and  45  to  rewrite  equations  41  and  42  as 

(vu  ■  N,  Vv  ■  N,  Vw  ■  Nj  ■  V  =  ^ a  (51) 

and 

(Vu  ■  N,  Vv  ■  N,  Vw  ■  ■  T2  =  ~j^P  (52) 

where 

a  =  (Vu  •  V,  Vv  ■  N,  Vw  •  JV)  •  fi  +  •  fl7Vv  •  f1;  Vw  -T^-N  (53) 

and 

P  =  (vu  ■  N,  Vv  ■  N,Vw  •  V)  •  f2  +  (Vu  •  f2,  Vv  •  f2,  Vw  ■  f2)  •  N  (54) 
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with  the  “haf  superscript  defined  as  outlined  above. 

Finally,  equations  44,  45,  47,  51,  and  52  can  be  compiled  to  obtain 


N 

Ti 

T2 


[V«] 

[Vu] 

[Vw] 


0  0  0  \ 

a  0  0 

/3  0  0  / 


or  more  simply 


/  [llX\  [Uy]  [llZ] 

N  M  [vz] 

\  K]  [Wy]  \wz] 


N  \  /  0  0  0 

d  0  0 

T2  )  \  p  0  0 


(55) 


(56) 


Alternatively,  equations  44,  45,  and  47,  can  be  compiled  to  obtain 


+ 


(57) 


or 


/  N  \  /  [/iVii]  \  /  N  \  /  N  \  /  Vu 

I  V,  J  I  [/iVu]  J  j  7j  J  =  [//]  J  7)  J  (  Vu 
\  t2  J  \  [dVw]  /  \  T2  J  \  T2  J  \  Vw 


+ 


TV 

0 

0 


Vu  \ 
Vu 

Vw  y 


_»  T  -* 

(  n  \  /  o 

o  -  [/-/]  7\ 

V  0  /  V  T2 


(  Vu 
Vu 
y  Vw 


(58) 


using  equations  41  and  42  as  well.  This  can  be  rewritten  as 

t  -> 

(ll-iux\  [ pLUy ]  [nuz\  \  /  Vu  \  /  o  \  /  0  \ 

[pvx]  [/ Wy]  [pvz]  I  =  [/j]  I  Vu  J  I  Tx  J  I  f,  I  + 

[mwx]  [/iw.y]  [pwj  /  \  Vw  /  \  f2  J  \  T2  / 

/  Vu  \  /  0  \  T  /  0  \  /  Vu  \  T 

[/i]  TVTiV  Vu  TVtTV  -  [/It]  fi  fi  Vu  NtN  (59) 

\  Vw  J  \  T2  )  \  T2  )  \  Vw  / 


noting  that  the  right  hand  side  of  this  equation  only  involves  derivatives 
that  are  continuous  across  the  interface  as  opposed  to  equation  56. 
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Figure  3:  Steady  State  Air  Bubble  -  Boundary  Condition  Capturing 
Method 


Figure  4:  Steady  State  Air  Bubble  -  Delta  Function  Method 
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In  viscous  flows,  the  velocity  is  continuous  across  the  interface  implying 
that  the  material  derivative  or  Lagrangian  acceleration  is  continuous  as 
well.  That  is, 


Du 

Dv 

Dw 

Dt 

Dt 

Dt 

are  valid  jump  conditions  allowing  one  to  write 


Px 

(2  pux)x  +  ( p{uy  +  vx))y  +  (p(uz  +  Wx))z 

.  p . 

P 

'Py 

(mK  +  Vx))x  +  (2  pVy)y  +  (p(vz  +  Wy))z 

.  P 

p 

(60) 


(61) 


(62) 


'Pz~ 

(p(uz  +  Wx))x  +  ( p{vz  +  Wy))y  +  (2pwz)z 

.  p . 

P 

based  on  equations  33. 

For  details  on  how  these  jump  conditions  are  all  utilized  in  a  sharp  inter¬ 
face  approach  to  multiphase  incompressible  incompressible  flow  with  vis¬ 
cosity  and  surface  tension,  we  refer  the  reader  to  [6].  Notably,  they  showed 
that  this  sharp  interface  approach  can  dramatically  alleviate  parasitic  cur¬ 
rents  because  the  pressure  profile  is  not  smeared  out  over  a  characteristic 
width.  In  a  test  example  of  a  stationary  air  bubble  in  water  (without  grav¬ 
ity),  they  showed  that  the  boundary  condition  capturing  sharp  interface 
approach  reduces  the  parasitic  currents  by  a  factor  of  1000.  Figure  3  shows 
the  sharp  pressure  profile  captured  by  the  boundary  condition  capturing 
approach,  as  opposed  to  the  smeared  out  delta  function  approach  shown 
in  figure  4. 


3.2  Treating  Viscosity  Implicitly 

In  [6],  the  explicit  treatment  of  the  viscosity  imposed  a  stringent  time  step 
restriction  of  At  ~  0( Ax2).  This  poses  a  severe  limitation  for  high  viscosi¬ 
ties,  and  limits  the  spatial  size  and  time  scale  of  the  problems  that  can  be 
considered.  A  first  step  towards  an  implicit  treatment  of  viscosity  for  sharp 
interface  methods  was  proposed  in  [15].  Assuming  p  and  p  are  spatially 
constant  within  each  phase  and  applying  V  •  V  =  0  allows  one  to  write  the 
viscous  term  in  33  as  V  •  (z/W)  where  v  =  p/ p.  Unfortunately,  equation 
59  couples  the  u,  v  and  w  terms  together  meaning  that  one  would  require 


14 


(a)  Delta  function 


(b)  Boundary  Condition  Capturing 

Figure  5:  Comparison  of  the  smeared  out  delta  function  approach  and  the 
boundary  condition  capturing  method  for  bubble  animation. 


an  implicit  method  that  applies  to  the  entire  coupled  system  of  equations, 
as  opposed  to  the  usual  component  by  component  approaches.  Since  [15] 
was  mainly  focused  on  computer  graphics  applications,  they  made  the  sim¬ 
plifying  assumption  that  all  the  viscous  fluxes  were  balanced  setting  the 
right  hand  side  of  equation  59  to  a  three  by  three  zero  matrix.  This  allowed 
them  to  fully  decouple  the  viscous  terms  into  three  separate  scalar  equa¬ 
tions.  The  new  scheme  was  applied  by  removing  the  viscous  terms  from 
equation  35,  renaming  the  result  of  equation  36  as  V** ,  and  integrating  the 
viscosity  implicitly  with  three  decoupled  equations  for  u,  v  and  w  of  the 
form 

u***  =  u **  +  A tV  •  {vVu***) .  (64) 

The  discontinuous  coefficients  v  are  handled  using  the  method  of  [4]  as 
outlined  in  section  2,  resulting  in  an  effective  viscosity  given  by  equation 
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15  of 


(65) 


u+9  +  v~{l  —  6) 

for  stencils  that  cross  the  interface.  Since  this  final  velocity  field  cannot  be 
expected  to  be  divergence  free,  [18]  advocated  once  again  using  equations 
37  and  36  to  make  the  results  of  the  three  implicit  solves  divergence  free. 

We  note  that  [21]  considered  a  related  problem  where  the  u,  v  and 
w  equations  are  coupled  together  in  the  case  of  spatially  varying  viscosity 
(i.e.  fully  spatially  varying,  not  just  piecewise  constant  across  an  interface). 
They  devised  a  semi-implicit  scheme  that  treated  the  terms  that  couple 
the  equations  together  in  an  explicit  fashion,  so  that  three  decoupled  scalar 
equations  could  be  considered  in  the  implicit  part  of  the  method.  A  similar 
approach  could  be  taken  here,  adding  the  jump  conditions  from  equation 
59  explicitly  so  that  three  decoupled  scalar  equations  could  be  considered 
in  the  implicit  part  of  the  method. 

Although  [15]  made  other  simplifying  assumptions  as  well,  such  as  ig¬ 
noring  the  jump  in  viscosity  in  equation  48,  they  carefully  considered  the 
jump  in  pressure  due  to  surface  tension  effects  in  that  same  equation.  In 
fact,  they  carefully  compared  the  boundary  condition  capturing  approach 
to  the  smeared  out  delta  function  approach  showing  dramatic  differences 
between  the  methods  as  shown  in  figure  5. 

4  Incompressible  Flame  Discontinuities 

The  boundary  condition  capturing  approach  was  extended  further  in  [17] 
where  a  numerical  method  was  developed  that  allowed  for  a  velocity  discon¬ 
tinuity  across  the  interface  as  well.  That  paper  focused  on  premixed  fuels, 
where  the  combustion  zone  separating  the  unreacted  fuel  from  the  reacted 
products  is  assumed  to  be  infinitely  thin,  allowing  for  the  process  to  be 
modeled  as  a  two  phase  inviscid  incompressible  flow  with  discontinuities 
in  velocity  and  material  properties  across  the  interface.  Unlike  previous 
methods  based  on  delta  function  formulations,  [17]  maintains  a  sharp  ve¬ 
locity  profile  across  the  interface  allowing  the  flow  to  be  fully  divergence 
free  in  each  subdomain.  In  contrast,  the  delta  function  approach  smears 
the  velocity  jump  across  the  interface  resulting  in  a  flow  with  compressible 
character  near  the  interface.  This  is  especially  important  since  the  inter¬ 
face  velocity  depends  on  the  local  velocity  of  the  unreacted  fluid,  which  is 
difficult  to  ascertain  when  the  velocity  is  nonphysically  smeared  out  to  be 
continuous. 

For  a  simple  contact  discontinuity,  the  interface  velocity  is  equal  to  the 
local  fluid  velocity,  i.e.  W  =  V.  Often,  only  the  normal  component  of  the 
interface  velocity  is  required,  i.e.  W  =  DN  where  D  =  Uy  =  V  ■  N.  For 
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reacting  flow,  the  interface  moves  at  the  local  velocity  of  the  unreacted 
fluid  plus  the  flame  speed,  S,  which  gives  the  rate  of  conversion  of  the  un¬ 
reacted  material  into  the  reacted  material.  This  accounts  for  the  movement 
of  material  across  the  interface.  If  we  denote  unreacted  and  reacted  mate¬ 
rial  properties  with  ”u”  and  ”r”  subscripts  respectively,  then  the  normal 
component  of  the  interface  velocity  is  given  by  D  =  ( Vn)u  +  S.  Note  that 
the  normal  component  of  the  velocity  is  discontinuous  across  the  interface. 
The  flame  speed,  S,  is  typically  defined  as  S  =  So  +  an,  where  So  and  a 
are  constants  characteristic  of  the  reaction  and  k  is  the  curvature. 

The  equations  for  conservation  of  mass  and  momentum  in  divergence 
form  are  given  by 

Pt  +  v  •  ( pV )  =  0  (66) 

(pV)t  +  X7-(pVV  +  pI)  =  0  (67) 

and  the  divergence  free  condition  is  V  •  V  =  0.  The  mass  and  momentum 
fluxes  through  an  interface  surface  element  (moving  in  the  normal  direc¬ 
tion  with  speed  D)  must  be  continuous  across  the  interface,  implying  the 
standard  Rankine-Hugoniot  jump  conditions  across  the  interface 

\p(VN  -D)}  =  0  (68) 


p(  VN  -  D)  (V  -  DN)  +  pN 


=  0. 


(69) 


If  T)  and  T2  are  the  local  unit  tangent  vectors  to  the  interface,  and  the 
mass  flux  in  the  moving  reference  frame  is  denoted  by 

M  =  Pr  ((Vjv)r  -D)  =  pu  {{VN)U  -  D )  (70) 

then  equation  68  can  be  rewritten  as  [M\  =  0,  and  equation  69  can  be 
written  in  terms  of  normal  and  tangential  components  as 

M  [VN]  +  [p]  =  0  (71) 


M  [VTl]  =  M  [Vt2]  =  0. 


(72) 


When  D  ^  Vn,  then  M  ^  0  and  equation  72  becomes  [VrJ  =  [Vr2]  =  0 
showing  that  the  tangential  components  of  the  velocity  must  be  continuous 
across  the  interface.  When  D  =  Vjv,  as  in  the  case  of  a  contact  discontinu¬ 
ity,  M  =  0,  and  the  tangential  velocities  are  completely  uncoupled  across 
the  interface.  The  jump  in  the  normal  velocity  is  derived  as 


0  =[D] 


pVn  —  p{Vn  —  D) 
P 


{ VN  ]  -  M 
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so  that 


■ 

'1' 

V 

=  M 

_P_ 

which  combined  with  equation  71  implies  that 


(73) 


(74) 


The  velocity  jump  conditions,  73,  are  applied  to  extend  the  velocity  field 
in  each  region  across  the  interface.  For  example,  the  unreacted  velocity  field 
is  extended  into  the  reacted  region  by  defining  unreacted  ghost  velocities 
as 


=  ur  —  M 


(75) 


=  vr~  M  ( - - n2  (76) 

\  Pr  Pu  J 

and 

wu  =  wr  ~  M  (  — - -  J  n3  (77) 

V  Pr  Pu  J 

where  n i,  n2,  n3  are  the  components  of  the  interface  normal  computed  at 
the  location  of  the  appropriate  velocity  component.  These  ghost  values 
are  used  in  any  discretization  of  the  unreacted  fluid  velocity  which  crosses 
the  interface.  Similarly,  reacted  ghost  velocities  are  computed  in  a  band 
on  the  unreacted  side  of  the  interface  by  adding  the  jump  condition  to  the 
local  unreacted  velocity,  and  are  used  in  the  discretization  of  the  reacted 
fluid  velocity.  This  avoids  combining  the  discontinuous  velocities  across 
the  interface.  For  example,  in  equation  35,  the  intermediate  velocity  field 
V*  is  computed  for  both  the  real  fluid  values  and  the  ghost  fluids  values, 
so  that  the  ghost  fluids  values  can  be  used  in  the  discretization  of  the  right 
hand  side  of  equation  37,  avoiding  combination  of  the  intermediate  reacted 
and  unreacted  velocities  across  the  interface.  Also,  when  solving  equation 
37  for  the  pressure,  the  jump  in  pressure  given  by  equation  74  is  treated 
using  the  boundary  condition  capturing  technique  outlined  in  section  2. 

Consider  two  flames  both  with  speed  5  =  1  initially  located  at  x  =  —  .5 
and  x  =  .5.  The  unreacted  material  is  at  rest  in  the  center  of  the  domain. 
Dirichlet,  p  =  0,  boundary  conditions  are  specified  at  both  ends  of  the 
domain.  Initially,  the  reacted  velocities  on  the  left  and  right  hand  sides  of 
the  domain  were  specified  as  u  =  —4  and  u  =  4  respectively.  Figure  6  shows 
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velocity  (t=0) 


velocity  (t=.25) 


Figure  6:  Merging  flames. 


the  computed  velocity,  and  illustrates  the  ability  of  [17]  to  treat  merging. 
After  merging,  the  domain  contains  a  single  phase  incompressible  fluid 
which  must  have  a  constant  velocity.  In  the  case  of  compressible  flow,  a 
finite  speed  of  propagation  rarefaction  wave  would  lower  the  velocity  to  the 
average  of  the  two  reacted  velocities  (zero  in  this  case).  For  incompressible 
flow,  the  “rarefaction  wave”  moves  at  infinite  speed  and  the  velocity  drops 
to  zero  in  one  time  step  as  shown  in  the  figure. 

Consider  two  circular  flames  with  reacted  material  inside  the  circles  and 
the  unreacted  material  outside.  The  flame  speed  is  given  by  S  =  1  +  .01k, 
and  Dirichlet,  p  =  0,  boundary  conditions  were  used  on  all  sides  of  the 
domain.  Figures  7(a)  and  7(b)  show  the  velocity  fields  at  different  points 
in  time,  before  and  after  the  flame  fronts  merge,  respectively.  Note  that 
the  topological  change  (merging)  requires  no  special  treatment. 

[16]  exploited  the  method  from  [17]  to  produce  visually  realistic  ani¬ 
mation  of  low  speed  combustion  processes.  By  modeling  the  expansion  of 
the  fuel  as  it  reacts  to  form  hot  gaseous  products,  they  were  able  to  pro¬ 
duce  high  quality  animations  of  visually  full  flames.  Although  the  level 
set  based  premixed  flame  model  suffices  for  modeling  the  blue  flame  core 
where  the  reaction  is  taking  place,  the  blackbody  radiation  emitted  by  the 
hot  gaseous  products  gives  the  fire  its  yellowish-orange  color.  Furthermore, 
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velocity  field  (t=.02) 


(a)  Velocity  field  before  merging, 
velocity  field  (t=.035) 


(b)  Velocity  field  after  merging. 

Figure  7:  Merging  circular  flames. 
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as  the  temperature  cools  to  the  point  where  the  blackbody  radiation  is  no 
longer  visible,  smoke  and  soot  are  apparent  in  the  flame.  Realistic  visual 
depiction  of  the  black  body  radiation  and  smoke  and  soot  require  that  tem¬ 
perature  and  (smoke)  density  of  each  fluid  element  that  has  crossed  over 
the  reaction  zone  be  tracked  with  a  passively  advected  reaction  coordinate 
variable. 


5  Multiphase  Interacting  Flows 

Recently,  [18]  presented  a  level  set  method  for  the  simulation  of  multiple 
(more  than  two)  fluid  regions  with  differing  viscosities,  densities,  and  vis¬ 
coelastic  properties.  The  ghost  fluid  method  and  the  boundary  condition 
capturing  methods  of  [4]  and  [6]  were  used  to  treat  the  discontinuous  in¬ 
terfaces,  and  the  methods  from  [17]  and  [16]  were  used  to  incorporate  the 
ability  for  one  material  to  be  converted  into  another,  n  distinct  fluid  re¬ 
gions  are  represented  using  n  level  set  functions,  resulting  in  a  vector-valued 
level  set  function  at  each  point  in  the  domain,  <j){x)  =  {(j> \(x), . . .  ,<f>n(x)). 
A  novel  projection  method  was  introduced  to  decode  the  vector  of  level 
set  values,  providing  a  ’’dictionary”  that  translates  between  them  and  the 
standard  single-valued  level  set  representation. 

The  multiple  level  set  projection  method  ensures  that  the  following 
properties  hold  at  each  point  in  the  domain:  (PI)  If  (f)j  is  the  smallest 
element  of  </>,  it  is  the  only  negative  element  and  its  magnitude  represents 
distance  to  the  interface.  This  property  implies  that  4>i  is  a  signed  distance 
function  in  region  i  for  all  i,  and  that  each  point  in  the  domain  is  assigned 
to  exactly  one  region.  (P2)  If  PI  holds  and  (j>k  is  the  second  smallest 
element,  <pk  =  —<t>j ■  This  implies  that  the  level  set  for  the  region  that 
the  point  is  closest  to  but  not  inside  is  a  signed  distance  function  as  well. 
These  properties  are  consistent  with  the  standard  level  set  methodology.  In 
particular,  a  standard  level  set  function  </>  can  be  regarded  as  (j>  =  {<f> i,  fa) 
where  4>\  =  cf>  and  <f>2  =  — </>.  This  gives  a  dictionary  that  translate  between 
c t>  and  </>,  by  determining  <frj  and  (f>k  at  any  point  and  treating  them  as 
(f) i  and  (j) 2-  The  projection  algorithm  of  [18]  is  simply  the  following:  at 
each  point  in  the  domain,  the  smallest  two  elements  of  (f>  are  determined, 
and  their  average  is  subtracted  from  each  element  of  (f> .  Subsequently,  all 
geometric  information  can  be  computed  from  the  level  set  which  is  negative 
at  each  point.  When  level  set  values  are  needed  between  grid  points,  </>  is 
interpolated  at  the  desired  location  and  the  projection  method  is  applied 
to  the  resulting  vector  on  the  fly. 

The  particle  level  set  method  of  [22]  was  extended  to  multiple  level 
sets  in  [18]  as  well.  Each  level  set  has  an  associated  set  of  particles  that 
are  seeded  near  the  boundary  of  its  interior  region.  After  evolving  each 
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Figure  8:  Viscous  letters  splash  into  a  pool  of  water,  then  change  into  low 
density  inviscid  fuel  bubbling  up  and  burning  when  they  hit  the  surface 
(350  x  200  x  350  grid,  10  phases). 
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( t>i  and  all  the  particles  independently  in  time,  each  <pi  is  combined  with 
the  particles  by  rebuilding  (f>~  using  region  i’s  particles  and  (f>f  using  the 
particles  of  all  the  other  regions  combined.  As  in  [22],  particles  are  used 
to  correct  the  level  sets  after  advection  and  after  reinitialization.  The 
projection  method  is  applied  after  each  particle  correction.  Note  that  the 
projection  method  has  the  nice  property  that  it  preserves  signed  distance 
(unlike  the  method  proposed  in  [23]),  so  that  it  can  be  applied  after  the 
reinitialization  step  without  harm. 

Figure  8  shows  an  example  with  ten  separate  fluid  regions.  Initially,  the 
letters  are  heavier  than  the  water,  highly  viscous,  and  sink  to  the  bottom 
almost  as  if  they  were  rigid  bodies.  The  simulation  parameters  are  then 
changed,  lowering  the  density,  setting  viscosity  to  zero,  and  adding  surface 
tension.  This  makes  the  letters  bubble  to  the  surface.  To  illustrate  phase 
change,  the  letters  are  also  set  to  be  reactive  with  air  so  that  they  catch 
on  fire  when  they  reach  the  surface. 

6  Conclusion 

In  this  paper,  we  have  shown  a  variety  of  applications  of  the  boundary 
condition  capturing  method  for  the  variable  coefficient  Poission  equation 
proposed  [4]  (motivated  in  part  by  the  ghost  fluid  method  of  [1])  in  the 
field  of  computational  physics  as  well  as  computer  graphics.  Future  work  is 
bound  to  be  focused  on  more  accurate  methods  for  a  fully  implicit  viscosity 
treatment,  alleviating  the  time  step  restrictions  imposed  by  surface  tension 
forces,  and  general  improvements  in  the  order  of  accuracy.  However,  the 
work  to  date  has  proven  that  sharp  interface  treatments  can  be  used  on 
fairly  complex  phenomena  for  a  variety  of  application  areas. 
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